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I. INTRODUCTION 

The response of a manybody system to an external pertur- 
bation, particularly its collective mode response, forms an im- 
portant method of analysis in condensed matter physics. Since 
the experimental realisation of dilute gas Bose-Einstein con- 
densation (BEC) there have been several collective mode ex- 
periments H] 121 m m |5l, in which perturbations of the con- 
fining potential were used to excite the system. Of partic- 
ular interest is the 1997 experiment of Jin et al. at JILA 121, 
which determined the excitation frequencies for the lowest en- 
ergy quadrupolar collective modes over a temperature range 
spanning the condensation transition. At low temperatures, 
where the system was mainly condensate, the results were ac- 
curately described by simple meanfield theory |]6]|2llll. How- 
ever, the behavior of the collective mode frequencies at higher 
temperature, where a significant thermal fraction was present, 
proved much more difficult to describe. Indeed, the descrip- 
tion of these experiments has become the de facto standard for 
testing finite temperature quantum field theories of BEC, and 
has been largely responsible for the development of gapless 
(^1 and second order |10| theories of the trapped Bose gas. 
To date, the only fully quantitative theoretical descriptions 
of these results have been provided by the Zaremba-Nikuni- 
Griffin (ZNG) formalism calculations of Jackson et al. lllTl 
in 2002 and the second order theory of Morgan et al. lfT2l in 
2003. 

In Fig. [TJa) and (b) we show a schematic representa- 
tion of the quadupolar modes excited in the JILA experi- 
ment, characterized by the projection, m, of their angular mo- 
mentum onto the z axis. A large body of theoretical work 
has been conducted on the subject of the JILA experiments 
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and is nicely summarized in a recent review by Proukakis and 
Jackson [,24J . The temperature dependence of the m = 2 
mode was accounted for by gapless Hartree-Fock Bogoli- 
ubov (HFB) theory calculations in 1998 [131, which included 
anomalous average and manybody effects in the system de- 
scription (also see Refs. Il25l l26l l27l l28l ). However, gap- 
less HFB failed to account for the rather sudden upward shift 
in the m = mode frequency observed in experiments at 
T w 0.65Tc. An explanation for the unexplained behavior 
of the m ~ mode was first provided by Stoof and cowork- 
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Figure 1: Schematic view of the quadrupolar modes of an oblate 
condensate, (a) The m — excitation where the axial and radial 
widths of the condensate oscillate out of phase and (b) the m = 2 
excitations and the widths of the condensate in the x and y directions 
oscillate out of phase. 



ers |fT9l |29l (also see 1301 ). who suggested that it arose from 
the coupling of in-phase and out-of-phase oscillations of the 
condensate and thermal cloud. This hypothesis suggested that 
an adequate theoretical description would require a dynamic 
treatment of both the condensate and noncondensate parts of 
the system. The first such formalism was the ZNG finite tem- 
perature theory |31, 32, 33 1 in which the system description 
takes the form of a Gross-Pitaevksii equation for the conden- 
sate, coupled to a Boltzmann equation for the nonconden- 
sate. Jackson and Zaremba 111] applied the ZNG theory to 
model the JILA experiment and found relatively good agree- 
ment with the experimental results. The following year, IVIor- 
gan et al. lfT2l reported the results of a second order theory 
that were also in good agreement with the experimental re- 
sults. That theory, the culmination of seven years of work 
by Burnett, Hutchinson, IVIorgan, Proukakis and coworkers 
lIll[l0l[I2[l3l[l3[l5l[El[l2l[IS|34'l, consistently included 
the dynamical interactions between the condensate and non- 
condensate atoms. 

In this paper we develop the projected Gross-Pitaevskii 
equation (PGPE) formalism to model the experiment of Jin 
et al. im. The PGPE method is a c-field technique ll35l appli- 
cable to the study of finite temperature degenerate Bose gases. 
It includes interactions between low energy modes of the gas 
non-perturbatively and is applicable in the critical region, e.g. 
see |l36l|32l[38l|29J. Indeed, PGPE predictions for the shifts 
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Figure 2: Schematic view of the c-field and the incoherent regions 
for a Bose gas in a harmonic trap potential, and the approximations 
we employ in our treatment of the collective mode dynamics. 



II. FORMALISM 

We take our system to be described by the second quantized 
Hamiltonian 

H= I dv [ii\v)(-^\l'^ + V{v,t)\^{v) 
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[/o*^(r)^^(r)*(r)*(r) 



(1) 



where ^'(r) is the quantum Bose field operator, and = 
Airh^a/m is the interaction strength, with a the s-wave scat- 
tering length. The trap potential is given as 



V{v,t)^Vo{v) + 6V{v,t), 



(2) 



where 



Fo(r) = -m {ujIx" + + uj^) 



(3) 



in critical temperature BOl are in good agreement with experi- 
mental measurements |4i |. While this formalism has success- 
fully predicted equilibrium properties for a degenerate Bose 
cloud, there have been no quantitative comparisons to dynam- 
ical experiments, so our comparison to the experiments of Jin 
et al. |i2j| forms an important test of this theory. 

A central feature of the PGPE approach is the formal di- 
vision of the system modes into the classical region C 1,42] . 
which is simulated using the PGPE, and an incoherent re- 
gion I, for which a classical field treatment is inappropriate 
(see Fig. [2]i. The C region dynamics are accounted for in the 
PGPE description, and the I region dynamics could be treated 
using, e.g., a Boltzmann description. As previous theoreti- 
cal work has has shown, the full dynamical treatment of the 
noncondensate is crucially important in providing a correct 
description of the JILA experiments. However, a full dynam- 
ical treatment of the I region is a rather complex addition to 
the theory that we do not consider here. Instead, we simply 
ignore the dynamics of the I region, with the justification that 
many of the noncondensate modes exist in the C region and 
their dynamical effect is included in the PGPE. We critically 
analyze this approximation by quantifying the dependence of 
equilibrium and dynamic properties on the energy cutoff, tcut, 
which sets the division between the C and I regions. 

The organisation of the paper is as follows. In Sec. [ll|we 
review the PGPE formalism for equilibrium properties of a 
trapped Bose gas, before outlining our extensions to the the- 
ory to model the JILA collective mode experiments. The re- 
sults of our calculations for the equilibrium states, collective 
mode frequencies, and damping rates are presented in Sec. Ill 
In that section we also consider the relative phase between the 
condensate and noncondensate modes, and the cutoff sensitiv- 
ity of our predictions, before we conclude in Sec. |IV] The data 
used to prepare our initial states is summarized in Appendix 

m 



is the static harmonic trapping potential, and 5V{Y,t) is a 
time-dependent perturbing potential we discuss further below. 



A. Experimental procedure 

The theory we develop here is relevant to the finite tempera- 
ture excitation experiment undertaken by Jin et al. in Ref. | 2J. 
In that experiment a degenerate *^Rb Bose gas was prepared 
in a magnetic trap with frequencies = ^x.y = 27r x 129Hz, 
LOz = In X 365Hz, and initial temperatures ranging from 
0.4-1 .4rc- The total number of atoms increased with tem- 
perature, varying from about 5 x 10'^ to 60 x 10'' atoms over 
that temperature range, with a condensate number of about 
6 X 10^ ± 2 X 10^ for T < 0.9Tc. 

Two different symmetries of perturbation were investigated 
in experiments, chosen to effectively couple to the lowest en- 
ergy m — Q and m = 2 collective modes. To excite the 
collective mode the trap was perturbed for 14 ms and then 
evolved in the static trap for a variable hold time before the 
cloud was released and imaged after expansion (see Fig. |3ja)). 
The condensate and non-condensate components were deter- 
mined using bimodal fits to the absorption image, and the 
widths of each component were extracted as a function of 
time. These results were analyzed to give excitation frequen- 
cies and damping rates for both components. 



1. Time dependent perturbation 

The perturbation used to drive the m = mode was a weak 
sinusoidal modulation of the radial trap frequency (see Fig. 
[3jb)). For the m = 2 mode the trap frequencies in the x and 
y directions were modulated sinusoidally with tt phase differ- 
ence (see Fig. [3]^c)). For calibration, the dipole mode was also 
measured by centre-of-mass excitation (see Fig. [3jd)). 
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Figure 3: Experimental time sequence and form of trap perturba- 
tions used to excite the Bose gas. (a) Overview of time sequence used 
in the experiment to excite and observe collective excitations in the 
system, (b)-(d) Schematic representations of the various trap pertur- 
bations used (see text). Contours of equipotential in the a;j/-plane are 
shown for the unperturbed (solid lines) and for the perturbed (dashed 
lines) traps, (b) Symmetric perturbation used to drive the m = 
mode, (c) The perturbation used to drive the m — 2 mode corre- 
sponds to a rotating ellipse, (d) Trap center displacement used to 
drive the dipole mode. 



In the our approach to modeling these collective excitations 
we explicitly simulate the perturbation procedure used in ex- 
periments. To do this we use a perturbation potential of the 
form 



+ w^y^ cos{ujpt)} 



5V{r, t) = -^A{t) {^Ix^ cos{ujpt + ^ 

(4) 

where ujp is the perturbation frequency, is a phase factor be- 
tween the X and y perturbation, and A{t) is the dimensionless 
time dependent amplitude of the perturbation (see Fig. [Sja)) 
of the square pulse form 



A{t) = 



Aq, 0<t< 14ms, 
0, otherwise, 



(5) 



with Ao = 0.015. The choice of = (0 tt) in Eq. Q 
corresponds to the perturbation used in experiment to excite 
the m = (m = 2) mode. In experiment uip was chosen "to 
match the frequency of the excitation being studied", with the 
motivation that this should cause the system to oscillate at its 
natural frequency. 

To drive the dipole oscillation, we use a perturbation poten- 
tial of the form 

6V{r, t) = ^mujl{(f sm'^{ujxt) - 2xdsin{uJa:t)}, (6) 

where d = 0.034/xm is the amplitude of the sinusoidal motion 
of the trap in the x direction. 



B. PGPE formalism 

We briefly outline the projected Gross-Pitaevskii equation 
(PGPE) formalism, which is developed in detail in Ref. 
The Bose field operator is split into two parts according to 



^'(r) = VcW+^AiW, 



(7) 



where ipc is the coherent region c-field and i/j-i is the inco- 
herent field operator (see [35|). These fields are defined as 
the low and high energy projections of the full quantum field 
operator, separated by the energy Ccut, as shown in Fig. |2] In 
our theory this cutoff is implemented in terms of the harmonic 
oscillator eigenstates {(p„(r)} of the time-independent single 
particle Hamiltonian 



2m 



+ Voir), 



(8) 



i.e. e„iy9„(r) = iJoVn(r), with e„ the respective eigenvalue. 
The fields are thus defined by 



nec 



(9) 
(10) 



where the a„ are Bose annihilation operators, the c„ are com- 
plex amplitudes, and the sets of quantum numbers defining the 
regions are 



C = {n: e„ < Ccut}, 
I = {n : e„ > Ecut}- 



(11) 
(12) 



The applicability of the PGPE approach to describing the fi- 
nite temperature gas relies on an appropriate choice for Ccut^ 
so that the modes at the cutoff have an average occupation 
of order unity. This choice means that the all the modes in 
C are appreciably occupied, justifying the classical field re- 
placement a„ c„. In contrast the I region contains many 
sparsely occupied modes that are particle-like and would be 
poorly described using a classical field approximation. Here 
we treat these modes using a meanfield approach. 



C. Equilibrium states 

In this subsection we review our procedure for calculating 
finite temperature equilibrium properties of a trapped Bose 
gas. The basic approach is to treat the C and I regions as in- 
dependent systems in thermal and diffusive equilibrium. We 
discuss the treatment of these regions separately below. Fur- 
ther details on this procedure are given in Sec. 3 of 1351 . 



1. PGPE treatment of C region 
The equation of motion for tpc is the PGPE 



ih- 



dt 



ifoV'c + Pc{C/o|V'c|Vc} 



(13) 
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where the projection operator 



rc{F{r)} = Vn{r) f dr v*„{r')F{r'), (14) 

nec 

formahses our basis set restriction of ipc to the C region. The 
main approximation used to arrive at the PGPE is to neglect 
dynamical couplings to the incoherent region |44|. 



An important feature of Eq. ( 13 i is that it is ergodic, so that 



the microstates ipc evolves through in time form a sample of 
the equilibrium microstates, and time-averaging can be used 
to obtain macroscopic equilibrium properties. Our basic pro- 
cedure for finding equilibrium states consists of evolving the 
PGPE with three adjustable parameters: (i) the cutoff energy, 
ecut, that defines the division between C and I, and hence the 
number of modes in the C region; (ii) the number of C region 
atoms, iVc; (iii) the total energy of the C region, Eq- The last 
two quantities, defined as 

Ec = f drrc(Ho+^\^Pc\^)^Pc, (15) 



Nc - / dv\M^)\\ (16) 



are important because they represent constants of motion of 
the PGPE ( [T3] l, and thus control the equilibrium state of the 
system. 

Given choices of these parameters, randomized initial states 
are constructed satisfying those constraints (see Sec. 3 of Ref. 
Il35l ). and are evolved according to Eq. ( [T3] l. Typically evo- 
lution times of order 10 trap periods are used for the system 
to relax towards equilibrium [45 ,1. Further evolution for times 
of order 10-100 trap periods are used to sample equilibrium 
microstates or to time-average equilibrium properties. 

To characterize the equilibrium state in the C region it is 
necessary to determine the average density, condensate frac- 
tion, temperature and chemical potential. Many of these quan- 
tities are also important for characterizing the I region (see 
Sec.|nC2]). 

The average density is obtained as a time average over the 
field microstates, i.e.. 



'^c(r) 



IV'c(r,- 



(17) 
(18) 
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where {tj} is a set of Mg times (after the system has been 
allowed to relax to equilibrium) at which the field is sampled. 
We typically use ~ 500 samples over 100 trap periods to per- 
form such averages. The notation nc(r) emphasizes this is 
the density contribution from atoms residing in the C region, 
with the obvious property that Nq = J drnc(r). For later 
convenience, we note that the momentum properties of the C 
region are easily evaluated using the momentum field (j>c{p), 
obtained by Fourier transforming the spatial field, i.e. 



(19) 



To find the condensate number, A'cond, in our equilibrium 
state, we use the Penrose Onsager definition |46|, that iVcond 
is given by the largest eigenvalue of the one-body density ma- 
trix 



GW(r,r') = (V£(r)Vc(r')), 



(20) 



which we are also able to evaluate as a time-average. 

Finally, using the Rugh method 1,47 J , we are able to deter- 
mine the temperature (T) and chemical potential (/i) by time 
averaging. We refer to Refs. ||36ll48l for additional details of 
this procedure. 



2. Meanfield treatment of I region 

The average properties of the incoherent region can be cal- 
culated from the one-particle Wigner distribution 



m(r,p) 

where 

eHF(r,p) = 



1 



(21) 



exp(/3[eHF(r,p) - /^]) - 1 ' 

Fo(r) + 2C/o(nc(r) + ni(r)),(22) 



P 

2m 



is the Hartree-Fock energy, and /i is the chemical potential. In 
this semiclassical description r and p are treated as continu- 
ous variables. However, care needs to be taken to ensure that 



Eq. (21 1 is only applied to the appropriate region of phase 
space spanned by the incoherent region, i.e. single-particle 
modes of energy exceeding ecut- In phase space this region 
is 



f^i = i r,p 



P 

2m 



+ Vo{v) > ecut 



(23) 



This allows us to calculate the incoherent region density 



ni(r) 



dpWi(r,p), 



and atom number 



iVi = j dvni{v). 



(24) 



(25) 



3. Obtaining equilibrium states 

The overall algorithm for generating equilibrium states 
within the PGPE formalism is summarized as a three step pro- 
cess: 

1. Using selected values of {Eq, Nq, £cut}, an appropriate 
randomized state is constructed and evolved according 
to the PGPE. Using time averaging Ncond,nQ{r),T, 
and /i are calculated. 

2. Using nc(r),r, and ^, the incoherent region is ana- 
lyzed, yielding rii(r), and hence the total atom number 



N = Nc + Ni. 



(26) 
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3. The values obtained for N, T, and Ncond are compared 
to the desired values and the values of {Eq, Nq, ecut} 
are adjusted before returning to step 1 . 

This process is quite time consuming since steps 1 and 2 can 
take of order a day to complete on commodity PC hardware. 
Finding desired initial configurations is made somewhat easier 
in the stochastic projected Gross-Pitaevskii equation (SPGPE) 
formalism Il35ll49ll50l which allows direct control of T and /i 
rather than Eq and Nq, although we have not used that here. 

Following this procedure, we have sampled equilibrium 
configurations with condensate occupation A^cond ~ 6000 ± 
2000 over the temperature range O.SlTc to 0.83Tc (see Sec. 



Ill A I. These samples of the equilibrium state are used as ini- 



tial conditions for the collective mode excitation procedure 
we discuss next. We examine attributes of the initial states, 
particularly the dependence on ecut, later in this paper. 



D. Formalism for dynamical modeling of collective mode 
excitation 



The fundamental approximation in our treatment of collec- 
tive mode excitation is to neglect the the dynamics of the I 
region, and their influence on the C region. In this approxi- 
mation the many-body dynamics is described by the PGPE 



= Ho^/jc+Vc{{SV{r,t) + C/ol^cP) ^c}(27) 



which differs from the PGPE used to generate equilibrium 
states by the inclusion of the perturbation potential. The moti- 
vation for considering only the C region is that many noncon- 
densate modes and their effect on the collective mode dynam- 
ics are included in C. We critically examine this approxima- 
tion later. 



1. Initial conditions 

We begin our simulations at t = when the perturbation 
potential is first applied (see Eqs. Q and (jSjl). The per- 
turbation potential, and the ensuring dynamics it generates, 
break the ergodicity of the PGPE for some period of time af- 
ter the perturbation has concluded (until the system rethermal- 
ized when the collective modes have damped out). Ensemble 
averages of the dynamical system thus need to be taken as an 
average over many trajectories. For each trajectory, we take 
as an initial condition 



4^'Hr,t-0)-^^'i(r,r,), 



Eq/ 



(28) 



where we have used the notation tp^'^ to represent the equilib- 
rium states generated for the time independent potential (i.e. 



the states i/'clr, Tj) appearing in Eq. ( 18 1), and tp^^ to rep 



resent the j-th trajectory for the simulation of the collective 
mode dynamics. The subsequent evolution of il^l^\r, t) is ac- 
cording to Eq. (27 1, which excites collective modes in the 
system . 



2. Observations and analysis 

In the analysis of the system dynamics we present in the 
next section we make extensive use of the line density, defined 
for the j-th trajectory as 



dydz\^[^\r,t)\' 



(29) 



This quantity, for a single trajectory, is itself of interest as 
the spatial integration corresponds to a spatial averaging of 
the system over the many modes in the C region, and there 
is some evidence that single trajectories of the PGPE can be 
compared to single experimental results. However, we will 
also be interested in the trajectory average calculated as 



(30) 



where we use Mt = 20 for the trajectory averaged results 
presented in this paper. 

We will also be interested in the momentum space equiva- 
lent line densities. 



(31) 
(32) 



We emphasize that these line densities only include contribu- 
tions from atoms in the C region. 

We can also use trajectory averaging to obtain other quan- 
tities, such as the coherent condensate component of the sys- 
tem. We do this by extending the one-body density matrix to 
the nonequilibrium case and evaluating it with trajectory av- 
eraging, i.e. 

Mt ^ 

GW(r,r',<) = ^(4^)(r,t)) i^'i\v\t). (33) 

Diagonalising G*^^^ (r, r', t) at each time we can obtain the in- 
stantaneous condensate (coherent) field ?/'cond(r, t), and hence 
the condensate line density 



(34) 



in. RESULTS 

In this section, we present a detailed analysis of the PGPE 
simulations of the JILA experiment |2| 



III A 



First, in Sec 

we present the parameters of the equilibrium states we have 
generated that we use as the basis for our collective excitation 
modeling. Then in Sec. |IIIB| we develop convenient observ- 
ables and examine the density response of the Bose cloud to 



the perturbative drive. In Sec. IIIC we present results for the 
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Figure 4: Equilibrium state properties, (a) Temperature, (b) total 
atom number, and (c) condensate number as a function of T/Tc. 
PGPE results (diamonds) and lines are guides to the eye. 



frequencies and decay rates of the lowest energy m — 2 and 
TO = modes. Then, in Sec. HID we analyse the effect of the 



energy cutoff in our formalism, and provide evidence for how 
it affects the equilibrium and dynamic properties of the Bose 
cloud. For completeness, we then calculate the frequencies of 



the dipole mode as a function of temperature in Sec. Ill E and 
lastly discuss the phase of the noncondensate and condensate 
oscillations in Sec. IIII Fl 



A. Equilibrium states 



First we present a summary of our results for the equilib- 
rium states generated according to the procedure discussed in 



Sec. II C 3 The macroscopic parameters of the states we have 
produced are shown in Fig. |4] These states provide initial 
conditions over the temperature range O.SlTc ^ O.SSTc with a 
condensate number in the range 3.5 x 10"^ — 7.5 x lO'^, which 
is comparable to the spread in condensate values used in ex- 
periment over this temperature range (see. Fig. l(c)ofRef. 
i2 |). A complete list of the parameters and properties of our 
initial equilibrium states is given in Appendix [A] 
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Figure 5: Line densities evolution after perturbation, (a) A single 
trajectory position space line density nf\x,t). (b) Trajectory av- 
eraged position space line density ni{x, t). (c) Condensate position 
space line density n^°"'*(a;, f). (d) A single trajectory momentum 
space line density rif\px-,t). (e) Trajectory averaged momentum 
space line density ni{px,t). Results for a system with T = 154 nK, 
iV = 2.0 X 10'', Ncond = 5.7 X 10^ 



B. Density response 

In this section, we show examples of the density response of 
the system after the sinusoidal perturbation has been switched 
off, and the cloud is evolving in situ in a static harmonic po- 
tential. 

Figures |5ja)-(e) and Figs. |6|a)-(b) show the evolution of 
the position and momentum line densities for a Bose gas af- 
ter the perturbation with to = symmetry has been applied, 
where the time is measured with t — corresponding to the 
beginning of the perturbation (see Fig. |3]l. The timescale of 
these results corresponds to the period of observation used in 
experiments. The position line density has a clear width os- 
cillation induced by the perturbation. We have made similar 
observations of the y and z line densities (defined analogously 
to Eqs. ( 29 1 and ( 3 1 1) and have verified that width oscillations 



also occur For the to = symmetry perturbation we find 
that the x and y oscillations are in phase, whereas the x and z 
oscillations are out of phase. Thus we conclude that the per- 
turbation has excited the to = mode more strongly than any 
other mode. 
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between the position space width and momentum space peak 
value for the condensate arises through the Heisenberg rela- 
tionship, i.e. the position and momentum widths of the con- 
densate mode are inversely related. The intermediate minor 
peak arises because of the out of phase oscillation of conden- 
sate width in the different directions integrated over to obtain 
the line density 1 54 1. 

In addition to the dominant condensate peak at — 0, a 
broad background feature is apparent in the momentum den- 
sity at larger \px\ values. This feature, which we attribute 
to the non-condensate portion on the system in the C region, 
is more clearly apparent in momentum line density shown in 
Figs. |7](a) and (b). 



C. Frequencies and decay rates of collective modes 

1. Observables 



Figure 6: Surface plots of trajectory averaged line density plots, (a) 
Trajectory averaged position space line density ni(x, t). (b) Trajec- 
tory averaged position space line density ni{px,t). Same data as 
displayed in Fig. |5[b) and (e) 
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Figure 7: (a) Single trajectory and (b) trajectory averaged momen- 
tum line density aX t — 32 ms for the condensate domain (dark 
shaded region) and the noncondensate domain (light shaded region) 
are shown (see text). Results for the same parameters given in Fig. 

m 



A similar study of the density response of the system to the 
perturbation with 171 — 2 symmetry reveals expected behav- 
ior: the widths in the x and y directions oscillate out of phase, 
and the z width remains (approximately) constant. 

Figure |5|d) shows the momentum line density for a single 
trajectory of a Bose gas after the perturbation with m = 
symmetry has been applied. Figure |5|d) and Fig. |6|b) show 
the trajectory averaged line density. The momentum line den- 
sity is sharply peaked at p^. = due to the presence of a con- 
densate. The peak value of the momentum line density oscil- 
lates periodically with minor peaks occurring between major 
peaks (see Fig. |6|b)). The major peak occurs first at t « 16 
ms and then returns each time the condensate width reaches 
the outer turning point of its oscillation in position space (i.e. 
the condensate is at its widest, see Fig. |5jc)). This connection 



It is necessary to measure appropriate observables to deter- 
mine the frequencies and damping rates of collective modes 
excited. In experiments bimodal fitting of the expanded sys- 
tem provided such observables, and gave independent infor- 
mation for the condensate and noncondensate (or thermal 
cloud). The in situ momentum distribution (e.g. see the mo- 
mentum line density in Fig. |7|b)) clearly reveals the distinct 
character of the condensate and noncondensate components, 
and (like in experiments) fitting a bimodal distribution to de- 
termine the widths of the condensate and thermal components 
would seem to be an obvious choice for observable. However, 
the in situ condensate momentum peak is extremely narrow 
and we have found that performing bimodal fits to the mo- 
mentum line density is ambiguous and noisy. We note that 
in experiments the expansion procedure gives rise to consid- 
erable broadening of the condensate momentum distribution 
(e.g. see ifSTi) and thus cannot be compared directly to our in 
situ line density. 

However, we can develop two useful observables that avoid 
the need for fitting. First we observe that the two momentum 
domains, defined as 



C = {Px ■■ \Px\ < Po}, 

n = {Px ■■ \Px\ > Po}, 



(35) 
(36) 



with Po = ^/2hmujx, are dominated by the condensate (i.e. 
narrow peak) and noncondensate (broad background) respec- 
tively (see shaded regions in Figs. |7] (a) and (b)). We thus 
refer to these domains as the condensate (i.e. c) and noncon- 
densate (i.e. n) domains respectively. The value of po is in 
some sense arbitrary as long as it is greater than the conden- 
sate momentum width, and much less that the characteristic 
thermal momentum (pth = h/X^h, with Adb the thermal de 
Broglie wavelength). Our choice, po = y/2hmujx, satisfies 
both of these these criteria. 

We can now define our two observables, as the variance of 
the momentum line densities on these restricted domains, i.e. 



Pcit) 



{pD 
(pI) 



(Px) 

(p.) 



(37) 
(38) 
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Figure 8: (a) Condensate domain observable, Pc{t), and (b) non- 
condensate domain observable, Pn{t), as a function of time. For 
a system excited by the m = symmetry perturbation: Trajec- 
tory averaged result (black line), and single trajectory result (grey 
line). For unperturbed equilibrium system: Trajectory averaged re- 
sult (black diamond line), and single trajectory result (grey diamond 
line). System parameters: T = 159 nK^ 0.75Tc, iV = 2.1 x 10*, 
iVcond = 4.9 X 10^ 



where 



(39) 



with cr = {c, n}. We note that the denominator of Eq. (39 1 
appropriately normalizes the moments, and the choice of vari- 
ance for Per, rather than the second moment (p^)(j, is to re- 
move the effects of any residual center of mass motion of the 
system. 

In Fig. [sj^a) and (b) we show examples of Pc(^) and Pn{t), 
evaluated from the PGPE simulation of an equilibrium system 
and a system excited by the perturbation with m — symme- 
try. From these results it is clear that the observables reveal the 
collective mode induced by the perturbation compared to the 
much smaller thermal fluctuations in the equilibrium states. In 
the collective mode analysis we always use Pc{t) and Pn{t) 
evaluated from the trajectory averaged line density, ni{px, t), 
however the results in Fig. [Slshow that if the single trajectory 
line density, n\^Hp.,t), is used to evaluate these quantities a 
useful signal is also obtained. 

For both Pc{t) and Pn(^) we notice that considerable damp- 
ing occurs over the period of observation. In both signals an- 
harmonic features are present, but are most apparent in the 
condensate observable where a weaker intermediate dip is ap- 



0.22 
> 0.2 

CM 

E 

E 0.18 



0.14 



0.12 



\ (a) A 
1 1 \^ l\ 

/ \ 1 \ A— A 














t [ms] 

Figure 9: Fits to observables: (a) Condensate domain observable, 
Pc{t), (line) and (b) noncondensate domain observable, Pc{t), (line) 
as a function of time for a system excited by the m = symme- 
try perturbation. Fits to observables using Eq. \40\ shown (dotted 
lines). System parameters: T = 159 nK^ 0.75Te, TV = 2.1 x 10'', 
Ascend =4.9 X 10^ 



parent. The origin of this feature is the same as for the in- 
termediate peak in Fig. |6]^b) (see discussion in Sec. IIIBi: 
Integration over the out of phase oscillation of the m = 
mode in the z direction. We have also verified that the observ- 
able signal is relatively insensitive to small adjustments of the 
value of pq used to define the c and n domains. 

As in the experiment we fit a decaying sinusoid of the form 



Pfit(i) = Ae-''* sm{ujt + i9) + B, 



(40) 



to our results, to obtain the collective mode frequency (uj) and 
damping rate (7). Example fits to the observable, shown in 
Figs. |9|a) and (b), indicate that while our combination of ob- 
servable and fitting function is adequate for accurately deter- 
mining the mode frequency, it does not provide a good de- 
scription of the amplitude or damping behavior of the modes. 



2. Mode frequencies 

Our results for the mode frequency variation with temper- 
ature are presented in Fig. 



10 along with the experimental 



results from Jin et al. \2 \ for comparison. We show results for 
the m = mode and the m = 2 mode, and give the frequen- 
cies for both the condensate (solid symbols) and nonconden- 
sate (open symbol) components. 
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Figure 10: Results for the frequency dependence on temperature for 
PGPE and experimental results of |2|. (a) m = mode frequencies. 
{h) m — 2 mode frequencies. Experimental results (grey symbols) 
and PGPE results (black symbols). Frequency for condensate (solid 
symbols) and noncondensate (open symbols). Error bars on some 
PGPE results indicate the spread in values from different calculations 
at the same temperature. 
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We first examine tlie m = mode beliaviour sliown in Fig. 
[TO^a). At temperatures below 0.6Tc the m — mode frequen- 
cies of the condensate and noncondensate components are al- 
most the same, indicating that the two components oscillate 
together In this temperature range the agreement with the ex- 
perimental results for the condensate frequency is good. There 
are no experimental measurements for the noncondensate be- 
havior in this regime as the noncondensate fraction is too 
small to measure. At temperatures above 0.6Tc our theoret- 
ical predictions and the experimental results exhibit markedly 
different behavior: As temperature increases above O.GTc our 
results (for both the condensate and noncondensate) decrease 
in frequency, whereas the experimental results show a rather 
rapid increase in frequency. This feature of the experimen- 
tal results evaded theoretical description (e.g. see |14|) until 
the works of Jackson et al. 1 1 1 1 in 2002 and Morgan et al. 
IIT2II in 2003. We discuss the origin of the disagreement be- 



tween PGPE and the experimental results further in Sec. Ill D 



and show that it arises from our lack of a dynamical descrip- 
tion of the I region. We note that the PGPE predictions of a 
downward trend in the frequency of the m — Q mode is con- 
sistent with the results of gapless Hartree-Fock-Bogoliubov 
calculations (see Fig. 2 of Hutchinson et al. llT3l ). indicating 
that anomalous average effects are included in our description. 
Our predictions are also in good agreement with the second 
order theory of Morgan et al. for the mode frequency in the 
absence of direct thermal driving (see diamond symbols on 
Fig. l(a)ofRef. |[12l). 
In Fig. 



10 b) the m ~ 2 mode is considered. Here we 



see reasonable agreement between the PGPE predictions for 
the condensate oscillation frequency and the experimentally 
measured values at all temperatures simulated. At high tem- 
peratures our predictions lie slightly above the experimentally 
measured values in a similar manner to the full second or- 
der predictions (i.e. including thermal driving) of Morgan 
et al. (see open circles in Fig. 1(b) of Ref. ifTSl ). For the 
TO = 2 mode noncondensate oscillation frequency, we see 
poor agreement with experimental results. There are no other 
theoretical predictions for the m = 2 thermal modes for us to 
compare against as neither Ref. ifTTTl or IIT2I present results for 
this case. The PGPE predictions for the noncondensate mode 
at temperatures above 0.70Tc show that this noncondensate 
decouples from the condensate, and that its frequency is well 
above that of the condensate. This behavior is qualitatively the 
same as that seen in experimental results (with experimental 
results only available at temperatures above 0.78Tc), however 
the upward shift of the thermal mode frequency we calculate 
is much lower than that observed in experiments. We discuss 
the origin of this quantitative disagreement between the PGPE 
and experimental results further in Sec. 



HID 



Figure 11: Results for the damping rates variation with tempera- 
ture for PGPE and experimental results. Results are given for the 
condensate m = modes (triangles) and m = 2 modes (circles). 
Experimental results (grey symbols) and PGPE results (black sym- 
bols). Error bars on some PGPE results indicate the spread in values 
from different calculations at the same temperature. 



3. Mode damping 



In Fig. 11 we present the PGPE predictions for the damp- 



ing rates of the to = and m — 2 condensate modes, which 
we compare against the experimental results. Although there 
is considerable scatter in the PGPE results, they appear to be 
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consistent with the experimental measurements. In particu- 
lar, we observe that in the temperature range O.STc — O.GTc 
the m = mode decays most rapidly (i.e. larger 7), while 
at higher temperatures the m = 2 mode gradually takes over 
with a larger damping rate, broadly consistent with the ex- 
perimental findings. We note that our choice of observable is 
more appropriate for determining mode frequency than decay 
due to the non-sinusoidal shape of the observable signal (see 
discussion in Sec. 



Ill CI 



This will lead to a systematic shift 
in our predictions for the damping rate and may be responsi- 
ble for the general downward shift of our results relative to the 
experimental measurements. In future work we will look into 
other observables to improve the accuracy with which we can 
analyse the mode damping rates. 



D. Cutoff dependence 

In this section we investigate the dependence of equilibrium 
and dynamic properties of the system on the energy cutoff 
(ccut) used in our simulations. 



1. Dependence of equilibrium states on 



To consider the effect of varying cutoff we follow the pro- 



cedure discussed in Sec. II C 3 to prepare an initial state with 
a cutoff of Ccut = and equilibrium parameters of 

T = 154nK= 0.74Tc, iVcond = 5.8 x 10^, and n^^in = 0.65 
where rimin is the mean occupation of the highest energy 
mode in the C region (i.e. the least occupied C region mode). 
The quantity nmin is an important indicator of the PGPE va- 
lidity, as it allows us to ensure that all the modes in C are 
appreciably occupied. 

To investigate the cutoff dependence we down-project the 
equilibrium microstates t/jQ of this system according to 



(41) 



where V' is the projector for the cutoff e^^^ < Ccut- The ef- 
fect of this projection is to reduce the size of C to a smaller 
region, C', and thus remove the occupation and energy of the 
modes lying between e^„^ and ecut- Since the constants of the 
motion, Ec (< Eq) and Nqi (< iVc) have changed it is inter- 
esting to investigate if the equilibrium properties of the down- 
projected state differ from the original state. To check this we 
evolve ipQi according to the PGPE (13 1 (on region C'), and 



analyze the thermal state that ?/'c' describes, after it is given 
time to thermalize. For the results we present here, we change 
e^ut from Abhujx to SOhujx. Over this range the number of C 
region modes decreases from 5706 (for the original state with 



Ecut = idhiUx) down to 1575 (for e^^ 



SOhcUx), i.e. the 



total number of C regions modes changes by a factor of 3.6 
between the cutoff extremes we consider 

In Fig. 12 we present results for the equilibrium properties 
of our down-projected states. Fig. [T2|a) shows the popula- 
tion (rimin) of the highest harmonic oscillator state present in 
the simulation as a function of the energy cutoff. As e^^^ is 
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Figure 12: Dependence of equilibrium variables of down-projected 
state on the energy cutoff e'^^t. (a) Average occupation of highest en- 
ergy C region mode rimin (diamonds), (b) temperature T (diamonds), 

(c) condensate number A'^cond (diamonds) and (d) total atom number 
TV (diamonds) and C region atom number Nc (squares). Dashed line 
in (a) is a linear fit to nmin in the variable l/e'^ut, solid lines in (b)- 

(d) are linear fits to the data. In (d) the shaded regions indicate the 
relative number of atoms in the C and I regions for each value of 



lowered, we see that the number occupying this highest state 
increases, in a manner consistent with the equipartition occu- 
pation of this mode (i.e. n^in ^ ^BT/e'^^^). 

Figures [T2|b)-(d) show the results for the macroscopic pa- 
rameters T, A'cond- ™d N, respectively, of the down- 
projected state. We can see that these paramenters (excluding 
Nq) do not vary systematically with e'^^^, and conclude that 
the equilibrium parameters of our PGPE simulations are not 
dependent on the energy cutoff. These are the first results we 
are aware of showing the insensitivity of classical field method 
predictions to cutoff Of course there are limits to how low we 
can take e^^j, since our C region must represent the conden- 
sate mode accurately which requires us to use a cutoff energy 
greater than the condensate chemical potential. 
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2. Collective mode dependence on Ecut 

Above we have shown that the equilibrium properties are 
insensitive to the cutoff defining the portion of the system in 
the C and I regions. In contrast we would expect that the 
PGPE theory for simulating collective modes, as developed in 
this paper, will show dependence on the cutoff. Fundamen- 
tally this is because the full dynamics of the C region are sim- 
ulated, while the population of the I region is neglected. Thus, 
in situations where the noncondensate dynamics are important 
the number of noncondensate modes included in C will have 
a direct effect on the dynamical observables of the system. 
It would therefore seem desirable to include as much of the 
noncondensate population in the C region as is possible, i.e. 
increase ecut- However, there is a limit to how high we can set 
Ecut- As discussed in Sec. II B we formally require that all 



the C modes are appreciably occupied for the classical field 
approximation to be a valid description of the Bose gas. For 
Ecut = 4:6hu!x we have n^in ~ 0.65, and so there is limited 
scope for using higher energy cutoffs. 

In the absence of a dynamical theory for the I region, ad- 
justing the value of Ccut allows us a mechanism by which to 
qualitatively investigate the role of the noncondensate dynam- 
ics in the collective mode dynamics. As we increase ecut from 
30hu!x to 46hujx the percentage of the total number of atoms 
in the C region increases from approximately 50% to 65% 
(seeFig.[T2ld)). 

We now investigate the frequency dependence of the ni — 
and m — 2 modes on the energy cutoff. Our procedure is the 




Figure 13: Frequency dependence of mode excitation on energy cut- 
off, for (a) m = mode and (b) m = 2. Condensate frequency 
(solid symbols) and noncondensate frequency (open symbols). Solid 
lines are linear fits to the condensate data, and dashed line in (a) is fit 
to the noncondensate data. 



same as in Sec. IIIC 2 except that we consider a single tem- 



perature of 0.74Tc and sample our initial conditions for the 
PGPE from the equilibrium states of varying e^ut (i-^- those 
used to average for the macroscopic parameters shown in Fig. 
12 1. From these simulations we determine frequencies of os- 
cillation of the condensate and noncondensate components. 



at a frequency of a; « 21.^2: (which the expected value for the 
noninteracting limit of a thermal cloud). In the intermediate 
with the results shown in Fig. [13]; a) and (b). The frequency of cutoff range, e'^^Jhuj^ E [38, 44], a combination of the con- 
densate dominated and noninteracting limit behaviors occur, 
and we were unable to fit a single frequency to these values. 
This suggests that there is a cutoff value, ~ AOfkOx, at which 
sufficient noncondensate is dynamically simulated for it to os- 
cillate independently of the condensate. 

We now make some observations, from comparison of the 
goodness of the PGPE description of the experimental fre- 
quencies in Fig. 10 to the cutoff analysis of the modes in this 
section. For the m = modes, both the condensate and non- 
condensate predictions are in poor agreement with experiment 
at 0.74Tc, and both are observed to be cutoff dependent [Fig. 
[T3ja)]. For the m — 2 mode we find that: (i) The condensate 
dynamics, which are in good agreement with experiment, are 



the noncondensate m = 2 mode for e'^^^jhiOx € [38,44] are 
omitted as a single frequency fit of sufficient quality cannot be 
found (see discussion below). 

The TO = modes show a dependence on e^ut' with the 
frequency of oscillation increasing as e^uj increases. The non- 
condensate frequency increases at a greater rate than the con- 
densate, which is consistent with the increase in the conden- 
sate frequency arising from it being driven by the noncon- 
densate component. Morgan has also seen this effect in his 
second order treatment by examining the influence of includ- 
ing thermal driving on the condensate mode (see Ref. y_8J). 
Our results clearly indicate that including the dynamics of all 
noncondensate atoms is crucial to obtain a condensate mode 
frequency that would be comparable with the experimental re- 
sults of Jin et at. 

The TO = 2 condensate mode shows no almost dependence 
on the energy cutoff, suggesting that the noncondensate com- 
ponent does not couple strongly to this motion of the conden- 
sate. However, the noncondensate to = 2 mode does show 
cutoff dependence: At low e^ut (small thermal component) the 
noncondensate oscillates at the same frequency as the conden- 
sate, while at high e^ut (l^ge thermal component) it oscillates 



independent of cutoff [Fig. 13 b)]. (ii) The noncondensate re- 
sults, which are in poor agreement with the experiment at high 
temperatures, are strongly cutoff dependent [Fig. [T3jb)]. In 
general, these observations lead us to expect that cutoff inde- 
pendent predictions of the dynamical PGPE theory are likely 
to be accurate in the absence of a dynamical theory of the I 
region, while cutoff dependent predictions are unreliable. In 
the latter case a dynamical theory of the I region is required. 
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Figure 14: Frequency of dipole mode as a function of temperature. 
Condensate (solid squares) and noncondensate (open squares). Error 
bars on some results indicate the spread in values from different cal- 
culations at the same temperature. Linear fit to the condensate (line) 
and noncondensate (dashed) results. 



E. Dipole Mode 

It is rigorously known that a harmonically trapped system 
will have a center-of-mass motion oscillation mode at the trap- 
ping frequency (Kohn mode) Ii52i . This mode is an important 
test of theory and was analyzed in experiment |2 | for the pur- 
poses of frequency calibration. Due to the presence of a pro- 
jector in the PGPE theory the Kohn mode is not a constant of 
motion (see Ref. |53 |) and so for completeness we investi- 
gate the dynamics of this mode here. To do this we use the 



PGPE (13 1 following the same procedure for setting up simu- 
lations as was done for the m ~ and m~2 modes, but with 
the dipole perturbation potential (|6|. To analyse our data, we 
study the first moments {px)a for a — {c, n} (see Sec. IIIC 1 



to provide observables for the condensate and noncondensate 



behavior. Our results, shown in Fig 14 indicate that the con- 
densate and noncondensate components both oscillate at ap- 
proximately l.Owr, as expected. 



F. Relative phase of condensate-noncondensate oscillations 

The relative phase of the condensate and noncondensate os- 
cillations has played a central part in the explanation of the 
sharp jump in the frequency spectrum of the m = mode. 
Stoof and coworkers fW/^Ol argued that the anomalous jump 
was caused by a transition from out-of-phase to in-phase os- 
cillations of the condensate and noncondensate components at 
high temperature. Morgan flSl lends support to this theory by 
calculating the relative phase between the oscillations of the 
condensate and noncondensate components, and shows that 
at moderate temperatures (^ O.STc) the components oscillate 
out of phase, whereas at high temperatures (~ 0.8Tc) they os- 
cillate in phase. This is consistent with the physical picture 
that a large noncondensate fraction oscillating at the nonin- 




Figure 15: Relative phase between the condensate and nonconden- 
sate collective modes for the different symmetries studied: m = 
mode (triangles) and m — 2 (circles). Error bars on some PGPE 
results indicate the spread in values from different calculations at the 
same temperature. Solid line is a linear fit to the m = data and 
dashed line is a linear fit to the m = 2 data. Same parameters as in 
Fig. [13] 



teracting frequency 2ujx couples strongly to the condensate 
TO = mode and drives it at this higher frequency. Morgan's 
results for the to = 2 mode show the relative phase between 
the components increases with increasing temperature up to 
about 0.85Tc, at which point a slight decrease is observed to 
begin. 

Here we follow the method of Morgan ifTSl closely. We cal- 
culate the phase difference between the two components using 
the first oscillation cycle after the perturbation is concluded. 
We find the relative phase by using the difference in minima 
of the two observable curves (see Fig. |9]l as a fraction of the 
half period of the condensate oscillation, to give a result from 

These results 



zero to TT. We present our results in Fig. 15 



are in reasonable qualitative agreement with those of Morgan 
at low temperatures, where the relative phases of each mode 
are increasing with temperature, with the to = 2 mode hav- 
ing a larger phase angle to the m — mode at any given 
temperature (c.f. Fig. 1 1 ifTSl ). However, generally our pre- 
dicted values for the relative phase are less than those calcu- 
lated by Morgan, and more importantly, we do not see the 
sudden reduction in phase angle for the to = mode as 
temperatures increases above ^ O.TTc- The likely explana- 
tion for our disagreement is that the fraction of nonconden- 
sate being dynamically simulated is not great enough so that: 
(i) At moderate temperatures (~ O.STc) the noncondensate 
is being dominated by the condensate oscillation, leading to 
a smaller than expected relative phase between the compo- 
nents (for both TO = and m = 2 modes), (ii) At higher 
temperatures (> 0.7Tc) the noncondensate component is in- 
sufficiently dominant to effectively drive the condensate back 
in-phase with its natural oscillation (applying only to the res- 
onantly coupled TO = mode). 



To further investigate these effects, in Fig. 16 we show the 
dependence of the relative phase on the energy cutoff e^^^. 
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Figure 16: Relative phase between the condensate and nonconden- 
sate collective modes the m = (triangles) and m = 2 modes (cir- 
cles) as a function of the down-projected energy cutoff e'^^^. Solid 
line is a linear fit to the m — data and dashed line is a linear fit to 
the m = 2 data. Same parameters as in Fig. [T3] 



and hence noncondensate fraction in the C region. In partic- 
ular, we consider a system at temperature 0.74Tc for various 
the same as was examined in Sec. 



(i.e. 



HID I. At this 



temperature Morgan predicts a relative phase of 0.257r for the 
m = mode (and that with increasing temperature this phase 
decreases) and a relative phase of 0.457r for the m = 2 mode 
(which remains approximately constant with increasing tem- 
perature, before starting to decrease at about 0.85Tc). Our re- 
sults for the dependence on e^^^ shows that the relative phase 
of both modes increase with increasing cutoff, although the 
m — mode does so more slowly than the m — 2 mode. The 
rn = 2 behavior indicates that as the noncondensate com- 
ponent being simulated increases (i.e. as e^^j increases) its 
phase, relative to the condensate, becomes more independent. 
We note that while the m ~ 2 condensate mode is cutoff in- 
sensitive, the relative phase between the m = 2 modes shows 
a dependence because the noncondensate mode does change 
character with e'^^^^ (see Sec. IIID2i. The slower rate of in- 
crease in the the relative phase of the m — modes with 



cutoff (as compared to the m = 2 modes, see Fig. 16 1 may be 



indicative of the resonant coupling between the components. 

At our maximum value of cutoff (Ccut — 4:5hujx), about 
50% of the noncondensate atoms are included in the PGPE de- 



scription, and we speculate that a complete dynamical repre- 
sentation of the noncondensate would lead to this mode driv- 
ing the condensate and the return to an in phase oscillation. 



IV. CONCLUSIONS 

We have presented a comprehensive study of the excitation 
spectrum of a Bose cloud at finite temperature by modelling 
the experiment of Jin et al. [2J with the PGPE formalism. Our 
results for mode frequencies are in good agreement with ex- 
periment and other theories up to about O.GSTc. At temper- 
atures above this our theory continues to provide a good de- 
scription of the m — 2 condensate mode. CuiTently our the- 
ory fails to predict the sudden increase in the frequency of the 
TO = condensate mode at temperatures above 0.65Tc. The 
origin of this failure in the current formalism is that we only 
provide a dynamical description for the portion of the noncon- 
densate in the C region. 

We have also examined the dependence of PGPE results on 
energy cutoff used to define the C region. Importantly, we 
demonstrated the insensitivity of the equilibrium predictions 
to energy cutoff. The study of cutoff dependence in the col- 
lective mode results clearly reveals the importance of the in- 
terplay between condensate and noncondensate components 
in the mode behavior, and suggests a new practical validity 
check for the PGPE theory: Dynamical predictions (in the ab- 
sence of a dynamical theory for the I region) should be veri- 
fied to be independent of the cutoff energy. 

The results of this first study with the PGPE give us great 
confidence that this theory is capable of providing a full de- 
scription of the JILA experiments. To do this would require us 
to implement a dynamical description of the I region, which 
we are currently pursuing. 
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Nc X 10^ 




A'^cond X 10^ 


N X 10'' 


T [nK] 


T/T, 




17 


4.64 


7.2 


3.95 


0.613 


70 


0.51 


0.89 


17 


4.54 


7.25 


3.74 


0.64 


78 


0.55 


1.12 


18 


4.54 


7.3 


3.74 


0.591 


71 


0.52 


0.94 


19 


4.54 


7.5 


3.67 


0.598 


74 


0.54 


0.89 


17 


5.05 


7.5 


4.24 


0.78 


90 


0.6 


1.26 


32 


11.9 


12.3 


7.46 


1.75 


127 


0.64 


0.89 


32 


11.4 


12.3 


7.05 


1.7 


126 


0.65 


0.84 


32 


11.1 


12.3 


6.8 


1.7 


126 


0.65 


0.79 


32 


11.9 


12.6 


6.95 


1.88 


134 


0.66 


0.83 


32 


10.9 


12.3 


6.26 


1.62 


126 


0.66 


0.81 


31 


10.6 


12.1 


6.31 


1.64 


128 


0.66 


0.86 


30 


10.6 


11.9 


6.52 


1.75 


133 


0.67 


0.91 


32 


11.9 


12.8 


6.83 


1.99 


141 


0.68 


0.84 


36 


11.9 


14.4 


5.6 


1.95 


143 


0.71 


0.81 


37 


11.9 


15 


5.16 


1.96 


147 


0.72 


0.7 


36 


11.9 


14.6 


5.64 


2.01 


149 


0.72 


0.85 


33 


12.4 


13.5 


7.12 


2.33 


157 


0.72 


1.02 


35 


12.4 


14.4 


6.22 


2.27 


158 


0.73 


0.91 


32 


12.4 


13.2 


7.41 


2.38 


160 


0.73 


1.09 


46 


13.1 


22 


5.77 


2.03 


154 


0.74 


0.67 


34 


12.4 


14.1 


6.57 


2.42 


163 


0.74 


0.99 


38 


11.9 


16 


4.92 


2.15 


159 


0.75 


0.88 


35 


12.1 


14.6 


5.72 


2.34 


163 


0.75 


1.02 


33 


12.4 


13.8 


6.89 


2.51 


168 


0.76 


1.01 


36 


11.9 


15.5 


5.26 


2.5 


174 


0.78 


1.05 


33 


12.4 


14.1 


5.98 


2.79 


179 


0.78 


1.15 


27 


10.1 


11.7 


6.32 


2.47 


172 


0.78 


1.44 


32 


12.4 


13.8 


7.06 


2.98 


186 


0.79 


1.29 


36 


11.9 


16 


4.67 


2.83 


186 


0.8 


1.06 


37 


11.9 


16.5 


4.62 


2.72 


185 


0.81 


1.03 


32 


12.4 


14.1 


6.86 


3.32 


197 


0.81 


1.32 


25 


9.59 


11.3 


6.32 


2.96 


190 


0.81 


1.79 


36 


11.9 


16.5 


5.14 


3.24 


202 


0.83 


1.14 



Table 1: A summary of the equilibrium state parameters used for the results reported in Fig. |4]and used as initial states for the results presented 
in Sec. |1I1 C| The first three columns give the parameters used to generate the initial states and the remaining columns give the macroscopic 
parameters determined for these states. 
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